home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-02 / wctunits.zip / MATH.PAS < prev    next >
Pascal/Delphi Source File  |  1991-07-30  |  7KB  |  318 lines
  1. unit math;
  2.  
  3. { Written by William C. Thompson - wct@po.cwru.edu }
  4.  
  5. { This unit was written to perform several basic mathematical calculations }
  6. { This unit automatically generates a random seed }
  7.  
  8. interface
  9.  
  10. const e=2.7182818284905;
  11.  
  12. function adjust(a:real):real;
  13. procedure quad(a,b,c: real; var x1,x2: real; var im:real);
  14. function relerror(observed,actual: real):real;
  15. function pow(a,b:real):real;
  16. function sign(r:real):integer;
  17. function rmax(r,s:real):real;
  18. function rmin(r,s:real):real;
  19. function imax(m,n:integer):integer;
  20. function imin(m,n: integer):integer;
  21. function log(x,base: real): real;
  22. function tan(a:real):real;
  23. function arcsin(x:real):real;
  24. function arccos(x:real):real;
  25. function arctan2(y,x:real):real;
  26. function gcd(m,n:longint):word;
  27. function lcm(m,n:integer):integer;
  28. function prime(n:longint):boolean;
  29. function rnd(a,b:real):real;
  30. function rnd2(a,b:real):real;
  31. function gaussian(mu,sigma:real):real;
  32. function distance(x1,y1,x2,y2:real):real;
  33. function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
  34. function hero(a,b,c:real):real;
  35. function factorial(n:word):extended;
  36. function stirling(n:real):extended;
  37. function combination(n,r:word):word;
  38. function permutation(n,r:word):word;
  39.  
  40. implementation
  41.  
  42. function adjust(a:real):real;
  43. { Adjusts angle to fit into (-pi,pi] }
  44. begin
  45.   repeat
  46.     if a<=-pi then a:=a+2*pi;
  47.     if a>pi then a:=a-2*pi
  48.   until (a>-pi) and (a<=pi);
  49.   adjust:=a
  50. end;
  51.  
  52. procedure quad(a,b,c:real; var x1,x2: real; var im:real);
  53. { Solves any quadratic equation.  If im=0, x1 and x2 are
  54.   two real solutions.  If im<>0, x1 and x2 are real parts
  55.   of two solutions and im is imaginary part.  The two
  56.   solutions would be x1 + im I and x2 - im I }
  57. var d,q: real;
  58. begin
  59.   im:=0.0;
  60.   if a=0 then begin
  61.     x1:=-c/b;
  62.     x2:=-c/b
  63.     end
  64.   else begin
  65.     b:=b/a;
  66.     c:=c/a;
  67.     a:=1;
  68.     d:=b*b-4*c;
  69.     q:=-b/2;
  70.     if d<0 then begin
  71.       x1:=q;
  72.       x2:=q;
  73.       im:=sqrt(-d)/2
  74.       end
  75.     else begin
  76.       if b<0 then q:=q+sqrt(d)/2
  77.       else q:=q-sqrt(d)/2;
  78.       x1:=q;
  79.       if q=0 then x2:=-b
  80.       else x2:=c/q;
  81.       if x2<x1 then begin
  82.         d:=x1;
  83.         x1:=x2;
  84.         x2:=d
  85.         end
  86.       end   { d>=0 }
  87.     end  { a<>0 }
  88. end;
  89.  
  90. function relerror(observed, actual: real):real;
  91. { Relative error }
  92. begin
  93.   if actual=0.0 then relerror:=abs(observed)
  94.   else relerror:=abs(observed/actual-1)
  95. end;
  96.  
  97. function pow(a,b: real):real;
  98. { Computes a^b }
  99. begin
  100.   if a=0 then
  101.     if b=0 then pow:=1                           { 0^0 = 1 }
  102.     else if b<0 then pow:=exp(b*ln(a))           { force error }
  103.     else pow:=0                                  { 0^x = 0 }
  104.   else if a<0 then
  105.     if abs(b)<1e-10 then pow:=1
  106.     else if relerror(b,round(b))<1e-8 then
  107.       pow:=(1-2*ord(odd(round(b))))*exp(b*ln(abs(a)))
  108.     else if (relerror(1/b,round(1/b))<1e-8) and odd(round(1/b)) then
  109.       pow:=-exp(b*ln(abs(a)))
  110.     else pow:=exp(b*ln(a))                       { force error }
  111.   else pow:=exp(b*ln(a))
  112. end;
  113.  
  114. function sign(r:real):integer;
  115. begin
  116.   if r<0 then sign:=-1
  117.   else if r>0 then sign:=1
  118.   else sign:=0
  119. end;
  120.  
  121. function rmax(r,s:real):real;
  122. begin  if s>=r then rmax:=s else rmax:=r  end;
  123.  
  124. function rmin(r,s:real):real;
  125. begin  if s<=r then rmin:=s else rmin:=r  end;
  126.  
  127. function imax(m,n:integer):integer;
  128. begin  if m>=n then imax:=m else imax:=n  end;
  129.  
  130. function imin(m,n:integer):integer;
  131. begin  if m<=n then imin:=m else imin:=n  end;
  132.  
  133. function log(x,base:real):real;
  134. { Computes log of x base base }
  135. begin
  136.   if (base=0) and (x=0) then log:=1
  137.   else if (base=0) and (x=1) then log:=0
  138.   else log:=ln(x)/ln(base)
  139. end;
  140.  
  141. function tan(a:real):real;
  142. begin
  143.   a:=adjust(a);
  144.   if abs(a-pi/2)<1e-10 then tan:=9.99e+37
  145.   else if abs(a-3*pi/2)<1e-10 then tan:=-9.99e+37
  146.   else tan:=sin(a)/cos(a)
  147. end;
  148.  
  149. function arcsin(x:real):real;
  150. begin
  151.   if x<0 then arcsin:=-arcsin(-x)
  152.   else
  153.     if x=1 then arcsin:=pi/2
  154.     else arcsin:=arctan(x/sqrt(1-x*x))
  155. end;
  156.  
  157. function arccos(x:real):real;
  158. begin
  159.   if x<0 then arccos:=pi-arccos(-x)
  160.   else
  161.     if x=0 then arccos:=pi/2
  162.     else arccos:=arctan(sqrt(1-x*x)/x)
  163. end;
  164.  
  165. function arctan2(y,x:real):real;
  166. { Computes angle of point (x,y) }
  167. var at2:real;
  168. begin
  169.   if x=0 then
  170.     if y>=0 then arctan2:=pi/2
  171.     else arctan2:=-pi/2
  172.   else begin
  173.     at2:=arctan(abs(y/x));
  174.     if x>0 then
  175.       if y>0 then arctan2:=at2
  176.       else arctan2:=adjust(-at2)
  177.     else
  178.       if y>0 then arctan2:=adjust(pi-at2)
  179.       else arctan2:=adjust(-pi+at2);
  180.     end
  181. end;
  182.  
  183. function gcd(m,n:longint):word;
  184. { Greatest Common Divisor }
  185. var k,l,r:integer;
  186. begin
  187.   m:=abs(m);  n:=abs(n);
  188.   if m=0 then gcd:=n else begin
  189.     k:=m;
  190.     l:=n;
  191.     while n<>0 do begin
  192.       r:=m mod n;
  193.       m:=n;
  194.       n:=r
  195.       end;
  196.     gcd:=m
  197.     end
  198. end;
  199.  
  200. function lcm(m,n:integer):integer;
  201. { Least common multiple }
  202. begin
  203.   if (m=0) and (n=0) then lcm:=0
  204.   else lcm:=abs(round((m*1.0)*n) div gcd(m,n))
  205. end;
  206.  
  207. function prime(n:longint):boolean;
  208. { Tests a number for primeness }
  209. var
  210.   p: boolean;
  211.   i, limit: word;
  212. begin
  213.   n:=abs(n);
  214.   if n in [0..1] then prime:=false
  215.   else if n in [2..3] then prime:=true
  216.   else begin
  217.     p:=false;
  218.     i:=2;
  219.     limit:=round(sqrt(n));
  220.     while (i<=limit) and not p do begin
  221.       if (n mod i=0) then p:=true;
  222.       i:=i+1
  223.       end;
  224.     prime:=not p
  225.     end
  226. end;
  227.  
  228. function rnd(a,b:real):real;
  229. { divides interval from a to b, inclusive, into 65535 values and
  230.   chooses one }
  231. begin
  232.   rnd:=a+random(65535)/65534*(b-a)
  233. end;
  234.  
  235. function rnd2(a,b: real):real;
  236. { Similar to rnd, but MANY more divisions (over 4 million) }
  237. begin
  238.   rnd2:=a+(65535.0*random(65535)+random(65535))/(sqr(65535.0)-1)*(b-a)
  239. end;
  240.  
  241. function gaussian(mu,sigma:real):real;
  242. { Produces a Gaussian distributed random variable with mean mu
  243.   and standard deviation sigma }
  244. var
  245.   r,v1,v2: real;
  246. begin
  247.   repeat
  248.     v1:=rnd(-1,1);
  249.     v2:=rnd(-1,1);
  250.     r:=sqr(v1)+sqr(v2)
  251.   until (r<1) and (r>0);
  252.   gaussian:=v1*sqrt(-2*ln(r)/r)*sigma+mu
  253. end;
  254.  
  255. function distance(x1,y1,x2,y2:real):real;
  256. begin
  257.   distance:=sqrt(sqr(x1-x2)+sqr(y1-y2))
  258. end;
  259.  
  260. function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
  261. { Finds equation for line y=ax+b between (x1,y1) and (x2,y2) }
  262. begin
  263.   if x2=x1 then begin
  264.     findline:=false;
  265.     exit
  266.     end;
  267.   findline:=true;
  268.   a:=(y2-y1)/(x2-x1);
  269.   b:=(x2*y1-x1*y2)/(x2-x1);
  270. end;
  271.  
  272. function hero(a,b,c:real):real;
  273. { Finds area of triangle with sides of length s }
  274. var
  275.   s: real;
  276. begin
  277.   s:=a/2+b/2+c/2;
  278.   hero:=sqrt(s*(s-a)*(s-b)*(s-c));
  279. end;
  280.  
  281. function factorial(n:word):extended;
  282. { computes factorial }
  283. var
  284.   i: word;
  285.   f: extended;
  286. begin
  287.   f:=1;
  288.   for i:=1 to n do f:=f*i;
  289.   factorial:=f
  290. end;
  291.  
  292. function stirling(n:real):extended;
  293. { computes factorial according to Stirling's approximation -
  294.   this may not be correct...  it should be, but ask your math professor }
  295. var
  296.   s: extended;
  297. begin
  298.   s:=1/n;
  299.   s:=s*exp(n*ln(n));
  300.   stirling:=s*exp(-n)/sqrt(2*pi);
  301. end;
  302.  
  303. function combination(n,r:word):word;
  304. { computes nCr }
  305. begin
  306.   combination:=round((factorial(n)/factorial(n-r))/factorial(r))
  307. end;
  308.  
  309. function permutation(n,r:word):word;
  310. { computes nPr }
  311. begin
  312.   permutation:=round(factorial(n)/factorial(n-r))
  313. end;
  314.  
  315. begin
  316.   randomize
  317. end.
  318.